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ABSTRACT 

The thermodynamic behaviour of self-gravitating A-body systems has been 
worked out by borrowing a standard method from Molecular Dynamics: the time aver- 
ages of suitable quantities are numerically computed along the dynamical trajectories 
to yield thermodynamic observables. The link between dynamics and thermodynam- 
ics is made in the microcanonical ensemble of statistical mechanics. The dynamics 
of self-gravitating A-body systems has been computed using two different kinds of 
regularization of the newtonian interaction: the usual softening and a truncation of 
the Fourier expansion series of the two-body potential. A particles of equal masses are 
constrained in a finite three dimensional volume. Through the computation of basic 
thermodynamic observables and of the equation of state in the P — V plane, new evi- 
dence is given of the existence of a second order phase transition from a homogeneous 
phase to a clustered phase. This corresponds to a crossover from a polytrope of index 
n = 3, i.e. p = KV~ 4 / 3 , to a perfect gas law p = KV -1 , as is shown by the isoenergetic 
curves on the P — V plane. The dynamical-microcanonical averages are compared to 
their corresponding canonical ensemble averages, obtained through standard Monte 
Carlo computations. A major disagreement is found, because the canonical ensem- 
ble seems to have completely lost any information about the phase transition. The 
microcanonical ensemble appears as the only reliable statistical framework to tackle 
self-gravitating systems. Finally, our results - obtained in a "microscopic" framework - 
are compared with some existing theoretical predictions - obtained in a "macroscopic" 
(thermodynamic) framework: qualitative and quantitative agreement is found, with an 
interesting exception. 

Key words: celestial mechanics, stellar dynamics; equation of state; methods: nu- 
merical 



1 INTRODUCTION 

Many astrophysical problems, ranging from cosmology down 
to the formation of planets, demand a knowledge of the phys- 
ical conditions for the existence of stable equilibria as well as 
for the appearance of instabilities in self-gravitating systems, 
either in the form of large gas clouds or of collections of con- 
centrated objects (JV-body systems). As a consequence, on 
this subject there is a long standing and important tradition 
of both theoretical and numerical investigations. In view of 
the aims of the present paper, there are three fundamen- 
tal reference works. The first one is an old and interesting 
paper (Bonnor, 1956) where the equation of state (P — V) 
was derived for a spherical mass of gas at uniform temper- 
ature in equilibrium under its own gravitation, and where 
the modified Boyle's perfect gas law was related with the 



onset of a gravitationally driven instability. The other two 
papers consider N gravitationally interacting particles in a 
spherical box which, under suitable physical conditions, can 
undergo a catastrophic state (the Gravothermal Catastro- 
phe) characterized by the absence of a maximum of the en- 
tropy (Antonov 1962) or, in a suitable range of parameters, 
can exist in a clustered phase characterized by a negative 
heat capacity (Lynden-Bell & Wood 1968). Already in this 
latter paper, then confirmed by subsequent work (Aaronson 
& Hansen 1972), it was surmised that - in analogy with or- 
dinary thermodynamics - self-gravitating systems undergo 
a phase transition. This hypothesis was also supported by 
the theoretical evidence of a condensation phenomenon in a 
system of self-gravitating hard spheres (Van Kampen 1964). 
Also numerical experiments on the A-body dynamics sup- 
ported this analogy through the evidence of the core-halo 
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formation (Aarseth 1963) or, in simulations of galaxy clus- 
tering (Aarseth 1979), through a phenomenology which is 
reminiscent of a phase transition occurring in atomic or 
molecular systems. More recent evidence of the existence of 
a phase transition in gravitational systems has been given 
for plane and spherical sheets (Reidl & Miller 1987; Miller & 
Youngkins 1998; Youngkins & Miller 2000) and for particles 
constrained on a ring (Sota, Iguchi & Morikawa 2000). 

Though the "negative specific heat paradox" - arisen 
by Lynden-Bell & Wood in their explanation of Antonov's 
gravothermal catastrophe - seems to faint the analogy with 
laboratory systems, Hertel and Thirring (Hertel & Thirring 
1971) showed that a negative specific heat - which is strictly 
forbidden in the canonical ensemble of statistical mechanics 

- can be legally found in the microcanonical ensemble as a 
consequence of the breakdown of ensemble equivalence; this 
ensemble inequivalence has been later confirmed for a sim- 
plified model of gravitationally interacting objects (Lynden- 
Bell & Lynden-Bell 1977)Q. 

All the results of the three above quoted papers were 
worked out in a macroscopic thermodynamic framework. 
Thus we can wonder if and how such a macroscopic phe- 
nomenology can be retrieved in the framework of a micro- 
scopic dynamical description. In other words, in the spirit 
of statistical mechanics, we can try to link the dynamics 
of the elementary constituents of a self-gravitating system 
with its large scale thermodynamics. (Throughout this pa- 
per we use "macroscopic" and "microscopic" in a relative 
sense, thus the microscopic level can be that of individual 
atoms or molecules of a gas cloud, as well as that of stars 
in a galaxy). This is actually one of the aims of the present 
paper, where we address the problem of the statistical me- 
chanical treatment of self-gravitating systems in order to 
understand whether the known thermodynamics can be de- 
rived from purely gravitational dynamics. The other aim of 
the present paper, closely related to the previous one, is to 
tackle the clustering instability, i.e. the breakdown of the 
homogeneous phase into a core-halo phase, in the attempt 
to describe it in analogy with the other - better understood 

- laboratory phase transitions. 

The paper is organized as follows. In Section £j we 
briefly discuss why a-priori a statistical mechanical treat- 
ment of gravitational iV-body systems runs into difficulties 
and, on the other hand, why these systems can be considered 
bona fide ergodic after suitable regularization; moreover, we 
present an unconventional regularization strategy of the ex- 
act Hamiltonian (transformed into a "Fourier-truncated" 
Hamiltonian) , leading to models that do not afford any 
practical advantage in the numerical treatment of the N- 
body dynamics but that are of great prospective interest 
for its analytic treatment in a statistical mechanical con- 
text and thus to apply concepts and methods that already 
proved powerful in the study of Hamiltonian dynamics. In 
Section pL after a quick glance at the integration method 
and parameters of the equations of motion and at the Mon- 
teCarlo simulations of the canonical ensemble averages, we 



* The inequivalence of statistical ensembles is not limited to grav- 
itational systems, but it occurs whenever the range of the poten- 
tial is comparable with the size of a system (Cipriani & Pcttini 
2000) and has recently attracted a lot of interest (Gross 1997). 



show how thermodynamics can be worked out through dy- 
namics with the aid of microcanonical statistical ensemble 
and by borrowing from Molecular Dynamics useful formu- 
lae already in the literature. We report the caloric curves 
(temperature vs energy), the specific heat and an order pa- 
rameter vs energy (suitably scaled with N), and a compar- 
ison among the outcomes of the simulations obtained with 
the Fourier-truncated Hamiltonian and those obtained with 
the standardly softened Hamiltonian. A somewhat familiar 
phenomenology in the numerical study of phase transitions 
is found mainly through the order parameter. A strikingly 
strong ensemble inequivalence is found: the canonical Mon- 
teCarlo averages do not keep even the slightest trace of the 
change of the specific heat from positive to negative and of 
the phase transition. Having confirmed the a-priori expected 
ensemble inequivalence, we discuss the reasons to consider 
the microcanonical ensemble as the good representative sta- 
tistical ensemble for self-gravitating iV-body systems and 
we present the dynamically worked out P — V equation of 
state. In Section W we compare our results with their theoret- 
ical counterparts obtained within a thermodynamic macro- 
scopic approach. Qualitative and quantitative agreement is 
found and everything fits into a coherent scenario, though 
with a remarkable difference between the clustering transi- 
tion and the gravothermal catastrophe. Finally, in Section 
|H| some conclusions are drawn. In Appendix A some basic 
definitions and concepts of the Gibbsian ensemble formula- 
tion of statistical mechanics are briefly recalled. In Appendix 
B a technical detail about the mean-field approximation is 
sketched. 



2 SELF-GRAVITATING N-BODY SYSTEMS 

In principle star clusters, galaxies, clusters of galaxies seem 
to naturally call for a statistical description of their dynami- 
cal behaviour. However, as the existence of negative specific 
heats reveals, there are some difficulties due to the very spe- 
cial nature of gravitational interaction. The gravitational 
interaction is always attractive, unscreened and of infinite 
range, therefore it is not stable, i.e. the potential energy 
U(ri, . . . , rjv) of N gravitationally interacting masses does 
not fulfil the condition U(ri, . . . , rjv) > —NB, with B a 
positive constant, whence the so-called lack of saturation: 
in the N — > oo limit the binding energy per particle and 
the free energy per particle diverge. This is also referred to 
as a breakdown of extensivity of these fundamental physi- 
cal quantities, at variance with what is familiar for ordinary 
laboratory systems. Moreover, from a rigorous mathemat- 
ical viewpoint, there is no equilibrium ground state: as N 
increases, at N/V = const, the energy per particle increases 
with N and cannot be defined in the thermodynamic limit 
The absence of an equilibrium state means that a gravita- 
tional system does not behave thermodynamically because 
standard thermodynamics does not apply to evolving sys- 
tems. However, a thermodynamic description is still possible 
for self-gravitating iV-body systems provided that they are 

t Though, from a physical point of view, the thermodynamic 
limit does not exist and has no meaning for self-gravitating 
systems. 
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not in a strongly unstable phase (Saslaw 1985), and in the 
case of slowly evolving system we can resort to statistical 
mechanics. 

Moreover, statistical mechanics gives a correct ther- 
modynamics if the thermodynamic potentials are extensive 
quantities. In a self-gravitating system only entropy in the 
microcanonical ensemble - at least on finite time scales and 
if the system is slowly evolving - is an extensive thermody- 
namic potentiaj^]; therefore, the results obtained in other en- 
sembles (grancanonical and canonical) are a-priori expected 
to be in disagreement with the results obtained in the mi- 
crocanonical ensemble. 

In what follows, systems of TV gravitationally interact- 
ing point masses will be considered that are described by 
the Hamiltonian function 

N N 

where Yi = (xi,yi,Zi). For the sake of simplicity we shall 
put G = 1 and mi — 1, i — 1, . . . , TV. We remark that 
this choice is convenient to keep the statistical and ther- 
modynamic properties unaltered (e.g. the non-extensivity) , 
at variance with what is implied by other possible choices 
(Heggie & Mathieu 1986). 

2.1 Dynamics and statistics 

Any physical phenomenon occurring in the system described 
by Hamiltonian (|l|) must have its origin in the dynamics. 
With the exception of those systems that can be treated with 
the methods of celestial mechanics (mainly perturbation the- 
ory), the dynamics of generic TV-body systems can be worked 
out only through numerical simulations. The numeric ap- 
proach, which is continually improved (Aarseth 1999; Mey- 
lan & Heggie 1997, and references quoted in these papers), 
provides the "experimental" counterpart of the theoretical 
investigation of these systems. 

In order to describe the physics of TV-body systems 
through a few relevant macroscopic observables, dynamics 
must fulfil the requirements of ergodicity and mixing to jus- 
tify the application of statistical mechanics. It is in general 
impossible for Hamiltonian systems of physical interest to 
ascertain whether ergodicity and mixing are rigorously ver- 
ified. However, a generic non-integrable and chaotic Hamil- 
tonian system with a large number TV of degrees of freedom 
can be considered ergodic and mixing, at least in a physi- 
cal sense. In fact, after a famous theorem due to Poincare 
and Fermi (Poincare 1892; Fermi 1923a, b), generic systems 
with three or more degrees of freedom are not integrable, 
i.e. there are no nontrivial invariants of motion besides to- 
tal energy. Only global invariants (like total momentum and 
angular momentum) , due to global symmetries of the Hamil- 
tonian (|l|) (as space translations and rotations), can exist. 
Thus, once the initial condition of a gravitational TV-body 
system is assigned, the representative point of the system 

t From the fundamental relation 1/T = dS/dE, being T propor- 
tional to the kinetic energy per particle and making use of the 
virial theorem, NT and E must have the same N dependence 
hence S is extensive. 



moves on a 6TV — 10 dimensional hypersurface of its 6TV- 
dimensional phase space. The lack of nontrivial integrals of 
motion (i.e. not related to global symmetries) entails that 
all the 6TV — 10 dimensional hypersurface of phase space is 
accessible to a trajectory issuing from any initial condition. 
A possible source of non-ergodicity would be the coexistence 
of regular and chaotic motions: observed for the first time in 
the Henon-Heiles model (Henon & Heiles 1964), it is in gen- 
eral implied by the Kolmogorov-Arnold-Moser (KAM) the- 
orem (Thirring 1978), which, however, has no practical rel- 
evance at large TV Thus, self-gravitating TV-body systems, 
after some suitable regularization to make finite the phase 
space volume, can be considered ergodic, so that time aver- 
ages of physically relevant observables can be replaced with 
suitable statistical averages computed on a given 6TV — 10 
dimensional hypersurface of phase space. 

The dynamical instability of the gravitational TV-body 
systems (Miller 1964) implied one among the first numerical 
evidences of the existence of deterministic chaos in Hamil- 
tonian dynamics. Since then, further numerical evidence of 
the existence of chaos in the self-gravitating systems has 
been provided by several authors (Kandrup 1990; Quinlan 
& Tremaine 1992; Goodman, Heggie & Hut 1993; Kandrup, 
Mahon & Smith 1994, and references therein; Cipriani & 
Pucacco 1994; Cerruti-Sola & Pettini 1995; Cipriani & Di 
Bari 1998; Miller 1999). On the other hand, chaotic dynam- 
ics in a many dimensional phase space implies a bona fide 
phase mixing (Casetti et al. 1999), which means that time 
averages of physical observables converge to their statistical 
counterparts in a finite time (whereas ergodicity implies an 
infinite time convergence). 

Therefore, the use of microcanonical statistical mechan- 
ics is naturally motivated by the reasonably good ergodic- 
ity and mixing properties of the dynamics of regularized 
self-gravitating TV-body systems (Cipriani & Pettini 2000). 
Taking into account the conservation of energy, of center of 
mass position and momentum, and of angular momentum, 
the microcanonical volume in phase space reads 

r N 

un{E) = / II dr ' d PiS(H(p, r) - E) (2) 

i— 1 

8^ Ti A p, - L)o< 3 > (£, P* - P ) 5(3) E ri - R) ' 

i i i 

In principle, by means of u>n(E), we can compute statisti- 
cal averages of any physical observable defined through a 
function A(p, r), and we can also compute the thermody- 
namics of a self-gravitating TV-body system, the link being 
made by the entropy defined as S = fcs logtjjv(-B); ks is the 
Boltzmann constant (see Appendix A). 

2.2 Regularized N-body Hamiltonians 

Let us now consider a system constrained in a finite volume. 
At variance with the customary choice of a spherical box 

§ KAM theorem requires extremely tiny deviations from integra- 
bility to imply the existence of sizeable regions of regular motions 
in phase space, moreover these deviations from integrability can- 
not exceed a threshold value which drops to zero exponentially 
fast with TV (Thirring 1978; Casetti et al. 1999) . 
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(Bonnor 1956; Antonov 1962; Lynden-Bell & Wood 1968), 
let us consider a cubic box of side length equal to L. The 
reason for such a choice is to explicitly break the rotational 
invariance of the Hamiltonian (|l|) and, in so doing, the mi- 
crocanonical volume lun(E) simplifies top] 



ujn{E) 



r N 

/ JJd ri d Pl 5(H(p,r)-E) 



(3) 



and thus we can borrow from the existing literature the an- 
alytic expressions, derived using lon(E) of Eq.(^), of some 
basic thermodynamic observables that are then used in our 
numerical computations. Working out anew the same ana- 
lytic expressions using u)n(E) of Eq.Q is a non trivial task 
beyond the aims of our present investigation. However, this 
is not a severe restriction if we refer to almost non-rotating 
systems whose angular momentum, even if conserved, is neg- 
ligible. The cubic box, allowing fluctuations of the total an- 
gular momentum, is thus equivalent to considering a whole 
ensemble of almost vanishing angular momentum systems. 

Even though the assumption of a confining box could 
seem unphysical, it is a simple way of idealizing different 
physical aspects which depend upon the chosen boundary 
conditions. In fact, the assumption of periodic boundary 
conditions is as if we took a fragment out of the bulk of 
a large system where - in the average - the number of par- 
ticles remains constant and small local density and energy 
fluctuations of the subsystem take place. In this case when 
a mass exits the box in a given direction, another mass with 
the same energy enters the volume from the opposite side to 
keep constant the energy and the number of particles. The 
assumption of reflecting boundary conditions amounts to 
mimic the presence of a halo of diffuse matter whose gravi- 
tational potential field would act to confine the system. Both 
assumptions about the geometric constraints of the system, 
besides the explicit breaking of rotational symmetry, also 
guarantee the finiteness of the configuration space volume 
over which the integration in Eq. Q is performed. 

In order to guarantee also the boundedness in momen- 
tum space, so that the whole integral in Eq.(|^) extends over 
a finite region of phase space, the two-body interaction po- 
tential must be regularized. 

We adopted two different kinds of regularization. The 
first one is the standard softening, adopted in numerical sim- 
ulations, with the replacement (e.g. see Binney & Tremaine 
1987) 



(4) 



+ if 



where 77 is a small softening parameter that bounds below 
the interaction potential and in so doing prevents the oc- 
currence of arbitrarily large values of the momenta. This 
regularization is local in space. 

The second regularization is nonlocal in space. It makes 
use of the Fourier representation of the Green function G(r — 
r') for the Poisson equation 



V 2 G(r - r') = -47n5(r - r') 



(5) 



^ In fact, both periodic conditions and reflecting boundary con- 
ditions destroy the conservation of P and R. 



in a cubic box of side L. In fact, one has G(r— r') = l/|r — r'| 
with the following Fourier development (Jackson 1975) 



r - r' 



32 \ - 

■kL 2.^ 



i(kix) sm(k m y) sin(fc„z) sin(kix') sin(k m y') sin(k n z') 



l,m,n—l 



(P + m 2 + • 



where ki — n l/L, k m = it m/L and k n = n n/L. Thus, 
Hamiltonian (^) can be exactly rewritten as 

N ri N 

2m" V 3 " + Vyl + P " ) ~ ~HT 2^ l 1 _ Oij)mimj 



E 



which, in an arbitrarily large volume, is completely equivP 
lent to ([[]). 

In order to regularize the two-body potential in Hamil- 
tonian (ll), one can truncate the Fourier expansion in Hamil- 
tonian (71) by summing I, m, n from 1 up to a finite number 

The two regularizations are a-priori inequivalent: the 
former pertains events (close encounters) localised in real 
space, whereas the latter makes use of collective coordinates 
(the Fourier modes) which are not localised in real space. 
By truncating the Fourier expansion in Hamiltonian (Q) we 
neglect all the dynamical details occurring at length scales 
smaller than the smallest spatial wavelength in the expan- 
sion. Loosely speaking, this is reminiscent of standard meth- 
ods in statistical mechanics, mainly in the context of the 
renormalisation group theory, where relevant and irrelevant 
degrees of freedom at a given length scale are Fourier modes 
with wavelengths above and below some cutoff respectively. 
In order to ascertain to what extent a truncated model still 
retains some physically relevant feature of the exact iV-body 
system (^), it is necessary to compare the outcomes of differ- 
ent truncations, i.e. different M m , and to make a comparison 
between the results obtained by simulating the dynamics as- 
sociated with Hamiltonian (jl|) and with a truncated version 
of the Hamiltonian (Q). 

There are different reasons for making numerically 
heavier the already heavy iV-body problem. First of all, the 
spatial coarse-graining considerably lessens the dramatic ef- 
fect that close encounters have on a reliable numerical com- 
putation of Lyapunov exponents (though they are not re- 
ported here), as well as on other observables of dynamical 
relevance; the frequency and the effect of close encounters 
in numerical simulations of iV-body systems are in this re- 
spect too much enhanced in comparison with the physical 
reality of practically collisionless systems as galaxies, clus- 
ter of galaxies and, to a lesser extent, star clusters. Second, 
the approach to equilibrium, as well as many other non- 
trivial dynamical properties, are by far better revealed by 
the use of collective coordinates, as a long experience in a 
different context has widely witnessed (Casetti et al. 1999). 
Third, truncated Hamiltonians out of (M) have the interest- 
ing property of naturally allowing a mean-field-like decou- 
pling of the degrees of freedom which is very interesting in 



(6) 



sm(kiXi) sm(kmlli) sin(fc n z 4 ) sin(kiXj) sm(k m yj) sin(fc n Zj) 
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view of analytical computations of both statistical mechan- 
ical properties and of chaotic properties of the dynamics in 
the framework of a geometric theory of Hamiltonian chaos 
(Casetti et al. 2000); moreover, such a mean-field- like repre- 
sentation of the Fourier truncated Hamiltonians leads to the 
definition of an order parameter that is useful to detect the 
occurrence of a phase transition. The regularization helps to 
smooth the local fluctuations on small spatial scales which 
do not significantly affect the macroscopic behaviour of the 
system, but cause a noisy variation of the averages on the 
time scales of interest. Finally, a wealth of models could be 
generated, even by retaining very few modes, which could 
reveal different aspects of the astonishingly rich dynamics of 
the exact iV-body problem. 

A generic truncated model Hamiltonian reads 

N 1 16 N 

1=1 j,j=l 



sm{kiXi) sm(k m yi) sin(k„Zi) sm(kiXj) sm(k m yj) sm(k n Zj) 



(l 2 + m 2 +n 2 ) 



£ 

/ ,m,n— 1 

(8) 

where we have chosen G = 1 and mi = 1, i = 1, . . . , N. 
From the associated Lagrangian function - being p X i = Xi, 
Pyi — Vii Pzi ~ Zi - the following equations of motion are 
derived 

32 \ -« io l,m,n 

Xi = V E 
l.m,n 

32 

v* = — 2-- 



32 " w lS (i) 

To S 7Tb—, — : T 2 7 — ttt cos{kiXi) s\a{k m yi) s\u{k n Zi) 
L z — ' (l 2 + m z + n 2 ) 

l ,m,n—l 
^2 TYlS^*^ 

■j? ^ y-j— — 9 ^ sm(kiXi)cos(k m yi)sm(k n Zi) 
■i 



.. _ 32 nSl 



/ 776 — i V~T on sin(kiXi) sin(k m yi) cos(k n Zi) 

— ' (l z + + n z ) 

L ,Tn,n—l 



(9) 



(i) 

i = 1,...,N, and where we have introduced S l = 
St tm ,n - sm(kiXi)sm(k m yi)sm(k n Zi), with Si lJn>n = 
X^jLi sin(kiXj) sin(k m yj) sin(k n Zj), to put in evidence that 
the same quantities Si,m tn enter all the equations of motion, 
what obviously simplifies the numerical computations. 



3 NUMERICAL RESULTS 

The phase space trajectories of an Hamiltonian system are 
constrained on a constant energy surface in phase space; 
therefore, time averages computed along the numerical so- 
lutions of the equations of motion, of both the Fourier- 
truncated system and of the softened version of the Hamil- 
tonian (Q), are generically expected to converge to micro- 
canonical ensemble averages^ (see Appendix A). Thus, in 
order to work out the thermodynamics of self-gravitating 
./V-body systems, we borrow from Molecular Dynamics the 



II Both regularizations do not qualitatively alter the chaoticity of 
the dynamics. 



formulae that link microscopic dynamics with macroscopic 
thermodynamics (Pearson, Halicioglu & Tiller 1985). 

At variance, the numerical estimate of canonical ensem- 
ble averages requires to construct suitable random marko- 
vian walks in the full phase space. Along such random tra- 
jectories (no longer constrained on any energy surface), the 
averages of physical observables converge to canonical en- 
semble averages provided that the recipe to generate the 
random walk is that of the so-called Metropolis importance 
sampling of the canonical ensemble weight. 

3.1 Numerical integration 

The numerical integration of the equations of motion (^|) and 
of the equations of motion derived from the Hamiltonian ([!]) 
with the replacement (Q) has been performed by means of a 
symplectic algorithm (McLahlan & Atela 1992). Some runs 
to check the reliability of the results have been performed 
using also a bilateral scheme (Casetti 1995). Symplectic inte- 
grators update the coordinates and momenta of an Hamilto- 
nian system through a canonical transformation of variables; 
for this reason, symplectic algorithms ensure a faithful rep- 
resentation of an Hamiltonian flow because, in addition to 
the conservation of phase space volumes and of the energy, 
they guarantee the conservation of all the Poincare inte- 
gral invariants of a system. Actually, there are interpolation 
theorems (Moser 1968; Benettin & Giorgilli 1994) stating 
that the numerical flows obtained through symplectic algo- 
rithms can be made as close as we may wish to the true 
flow of a given Hamiltonian. Though locally other integra- 
tion schemes can give more precise results in presence of 
close encounters by using different time steps for individual 
particles (Aarseth 1985), the non-symplectic nature of these 
algorithms might a-priori somewhat alter the frequency with 
which different regions of phase space have to be visited by 
an ergodic dynamics, at least for very long runs. On the 
other hand, we are just interested in long runs, so that time 
averages of the chosen observables display a good stabiliza- 
tion, and, in order to safely replace microcanonical ensemble 
averages with time averages, dynamics has to properly sam- 
ple the phase space. 

The dynamics of the direct iV-body system has been 
numerically computed using integration time steps At rang- 
ing in the interval (5 • 10 -5 — 10~ 3 ): the relative variations 
AE/E of the total energy were in the range (10 -10 ,10 -4 ) 
on integration runs of (10 6 , 10 7 ) time steps. The softening 
parameter r\ has been set to 0.01 and scaled as r\ — r/d, where 
d = D/N 1/3 , with D = min(N 2 /E, L). 

For what concerns the Fourier-truncated model, compu- 
tations have been done with J\f w — 5, 7, 10, i.e. a total num- 
ber of modes equal to 125, 343 and 1000 respectively. The 
initial conditions on the particle coordinates (xi,yi, Zi) have 
been chosen at random with a uniform distribution in the 
interval (0, L). The initial momenta have been chosen at ran- 
dom with a gaussian distibution of zero mean and a suitable 
variance to approximately match the desired initial value of 
the total energy. An opportune velocity rescaling allowed to 
precisely fix the initial value of the total energy. With inte- 
gration time steps At ranging in the interval (0.01,0.001), 
we got relative variations AE/E of the total energy in the 
range (10 -8 , 10 -6 ) with zero mean, i.e. without any drift, 
and on long integration runs of 10 6 time steps. For E < 0, 
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the system has been then let evolve for about 3tn, where 
the dynamical time is defined as to <x (N/\E/N\ 3 ^ 2 ) (this 
makes the system virialize); for E > 0, we let the system 
evolve until it attains a stationary state of dynamical equi- 
librium between kinetic and potential energies. 

By varying the side L of the box, the volume V has 
been varied. With both periodic and reflecting boundary 
conditions, the initial velocity of the center of mass has been 
set equal to zero. The initial total angular momentum can 
be made very small, but then it fluctuates because of the 
explicit breaking of rotational symmetry. The larger N, the 
smaller such fluctuations can be. 



3.1.1 The N-scaling of the results 

The number N of interacting bodies has been varied in the 
range 25 — 500, with some checks up to N — 2000. In the 
numerical simulations of extensive Hamiltonian systems, i.e. 
with short-range interactions so that energy and other phys- 
ical observables are additive, the comparison among the re- 
sults obtained by varying N is naturally made through den- 
sities: the values of the observables divided by N. 

In what follows, we shall mainly vary the energy E by 
keeping the ratio g — N/V constant. From a dimensional 
point of view, the potential energy U = — y^'j-"' l/ r y i g 
[U] = N 2 /L, and as L = V 1/3 = {N/g) 1/3 (g is constant), 
we get U ~ g 1/3 N 5/3 . As we shall see in the following, this 
actually suggests the correct way of scaling the results ob- 
tained for different values of N. 



3.2 MonteCarlo computations 

For standard Hamiltonians H = yj. \p\ + U(q), the weight 
e~® H splits into a multidimensional gaussian, originated by 
the kinetic energy term, and into a configurational part 
e -pu(q). p j s t ne inverse of the temperature. Only this lat- 
ter part is dealt with by the MonteCarlo algorithm. The 
system under consideration has to be at equilibrium so that 
the transition probabilities W(l — > 2) and W(2 —> 1) be- 
tween any two configurations "1" and "2" satisfy the con- 
dition W(l -> 2) = W(2 -> 1), the so-called detailed bal- 
ance. A MonteCarlo move consists of a random update of 
the coordinates. If the system lowers its energy with the 
proposed configuration update, then the new configuration 
is accepted, otherwise the transition probability V^({gi} — * 
{q[}) = Nexp{-/3{U(q' i ) - Ufa)]}, where M is the nor- 
malization factor, is compared with a random number £ 
uniformly distributed in the interval [0,1]; if W > C, then 
the new configuration is accepted, otherwise the old config- 
uration is counted once more. This is the essential of the 
Metropolis importance sampling algorithm (Binder 1979). 
The average acceptance rate is usually kept not far from 50 
per cent by adjusting the mean variation of the coordinates 
at each update proposal. The MonteCarlo estimate of the 
canonical average of an observable A is then simply given 
by 



Nmc ^ 

{<?;} 



3.3 Dynamical analysis of thermodynamical 
observables 

In deriving the thermodynamics of self-gravitating iV-body 
systems, one of our aims is to address the problem of the 
existence and of the characterization of the clustering phase 
transition. It is not out of place to remind here that the ba- 
sic thermodynamical phenomenology of a phase transition is 
characterized by the sudden qualitative change of a macro- 
scopic system when its temperature varies across a critical 
value T c . This sudden change is qualitatively due to collec- 
tive microscopic behaviours and is quantitatively reflected 
by the singular!**] temperature dependence of the most rel- 
evant thermodynamical observables. Changes of state, like 
melting, condensation and so on are examples of first order 
phase transitions, the spontaneous magnetization of a para- 
magnetic material when temperature is lowered below the 
so-called Curie temperature is an example of a second order 
phase transition^]. The usual framework of both theoretical 
and numerical investigations of phase transitions is mainly 
that of Gibbs' canonical ensemble (see Appendix A). 

Recently, the question of whether microscopic Hamil- 
tonian dynamics displays some relevant change at a phase 
transition has been addressed by several works (Antoni 
& Ruffo 1995; Caiani et al. 1997; Dellago & Posch 1997, 
and earlier works therein cited; Gross 1997; Casetti et al. 
2000; Cerruti-Sola, Clementi & Pettini 2000). The dynam- 
ical approach has its natural statistical mechanical coun- 
terpart in the microcanonical ensemble (see Appendix A). 
After a relaxation time which depends on the initial con- 
dition, the time averages of observables computed along 
numerical phase space trajectories provide good estimates 
of their microcanonical counterparts. In the following Sec- 
tions we adopt this dynamical approach to investigate if 
also in self-gravitating iV-body systems a phase transition is 
present and if the dynamics shows a corresponding qualita- 
tive change. 



3.3.1 Caloric curve 

The first goal of our numerical computations was to obtain 
the caloric curves of different self-gravitating systems re- 
sulting from different choices of the parameters. The caloric 
curve T(E) (temperature vs energy) represents a basic link 
between the microcanonical statistical ensemble and ther- 
modynamics, thus between dynamics and thermodynamics. 
From the entropy definition S(E,V, N) — fcs \oguiN(E) in 
the microcanonical ensemble, ujv being given by Eq.(^), 
temperature is derived as 



dS(E) 
dE 



(10) 



the sum being carried over the Nmc accepted configurations 



T(E) 

which, for systems described by standard Hamiltonians H = 
y\ ipf + U(q), yields the expression 

True singularities are possible only in the N — > oo limit. Here 
"singular" means that the larger N the sharper is the almost 
singular pattern of an observable as a function of temperature, 
tt According to the Ehrenfest's definitions, the first or the second 
derivative of the Hclmoltz free energy with respect to temperature 
is discontinuous at a first order transition point or at a second 
order transition point, respectively. 
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l) (K- 1 ), 



(11) 



where {K~ 1 ) u is the microcanonical average of the inverse 
of the kinetic energy K. The replacement of this average by 
means of the time average {K~ 1 )t provides the dynamical 
estimate of the temperature. 

With the equivalent definition of microcanonical en- 
tropyas S(E,V,N) = Icb log Qn(E) (see Appendix A), 
Eq.(U(« yields the more familiar expression 



(12) 



relating temperature with the mean kinetic energy per de- 
gree of freedom. In our numerical computations, the two 
temperatures have been always found almost coincident, well 
within the theoretically expected deviation of 0(1/N). Any- 
way, in the results reported here, T is computed according 
toEq.®. 

The Fourier-truncated model, due to the finite number 
of modes M w considered, underestimates the potential en- 
ergy with respect to the direct model. However, since an 
additive constant in the Hamiltonian (^) does not affect the 
dynamics, we can shift a-posteriori the energy scale by com- 
puting the average energy difference of a set of random con- 
figurations whose potential energy is computed according 
to both Eqs.(|) and (|). Thus, in the following, in order 
to compare by superposition the results of the two models, 
the values of the energy densities obtained from the Fourier 
truncated model have been shifted towards the energy den- 
sities of the direct simulation. 

In Fig. |l| the caloric curve is reported for the Fourier- 
truncated model (|]) with M w = 7, N = 50 interacting ob- 
jects and L = 50, so that g = 4 • 10~ 4 . The results obtained 
with periodic and reflecting boundary conditions are synop- 
tically displayed together with the outcomes of MonteCarlo 
simulations in the case of reflecting boundary conditions. In 
the case of MonteCarlo simulations, temperature is the in- 
put parameter and the average total energy is the outcome; 
anyway, since there is a one-to-one correspondence between 
temperature and the average total energy, the results can 
be reported without ambiguity also on the T vs E plane. 
The dynamical (microcanonical) caloric curves display the 
same qualitative features and are even quantitatively not 
very different one from another. The low energy branch has 
a negative slope, meaning that the specific heat is negative. 
At some value of the energy per degree of freedom, the slope 
of T(E) becomes positive and the curve bends towards an 
asymptotic straight line of slope 2/3, proper to a gas of in- 
dependent particles. 

Fig. ^ shows the effect of different truncations in Eq. 
(^): qualitatively, the results are satisfactorily stable, quan- 
titatively, the larger J\f w , the closer is dT/de (with e = E/N) 
to the newtonian slope —2/3. 

Notice that Figs. ^ and ^ display also a slope —2/3 of 
the negative-cy branches, in agreement with the prediction 
of the virial theorem, in the regime in which the effect of the 
box is negligible (large negative energy values). 

Fig. ^ and Fig. ^ summarize a lot of information. They 
show the caloric curves obtained for different values of N, 
for the newtonian TV-body simulations and for the Fourier- 
truncated model computed with the canonical MonteCarlo 
algorithm and with the microcanonical simulations. The 




E/N 



Figure 1. Caloric curve T(E). q = 4- 1CT 4 is the density N/L 3 . 
Here N = 50 and J\f w = 7. Starred squares refer to periodic 
boundary conditions, full circles refer to reflecting boundary con- 
ditions (box). Open circles represent MonteCarlo simulation re- 
sults for the Fourier-truncated model in the canonical ensemble. 



data collapsingrj is obtained by rescaling the energy with 
iV 5 / 3 and temperature with N 2 ^ 3 (notice that temperature 
is already an en ergy divided by N). The simple argument 
given in Section 3.1.1 appears correct. 

There is a very good agreement between the results 
obtained with the Fourier-truncated model and those ob- 
tained with the direct A r -body simulations. This is a very 
interesting result because of the relatively small number 
of Fourier modes considered, and in view of adopting the 
Fourier-truncated model for the analytic computation of sta- 
tistical mechanical averages. 

It is remarkable that the change of sign of the slope of 
the caloric curve implies the existence of a phase transition. 
In fact, T(E) can only change the sign of its slope either 
through a "V" -shaped pattern or through an "U" -shaped 
pattern: in the former case, the derivative dT/dE is discon- 
tinuous at the cuspy point so that the specific heat makes 
a discontinuous jump (first-order phase transition), whereas 
in the latter case the same derivative vanishes at the mini- 
mum so that the specific heat diverges (second order phase 
transition) . 

The comparison among the dynamical (microcanonical) 
and the MonteCarlo (canonical) caloric curves shows a re- 
markable fact: the canonical caloric curve seems indistin- 
guishable from that of a perfect gas, no reminiscence in the 
canonical results seems to be left of the feature shown by 
the microcanonical curves. As already mentioned through- 
out this paper, a negative specific heat is strictly forbidden in 

it This is standard jargon in the numerical study of phase tran- 
sitions, and it means that when different sets of results are com- 
pared after a suitable rescaling with the parameter that varies 
from a set to another, then all the results crowd on the same 
curve. 
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Figure 2. Fourier-truncated model. Comparison between the 
caloric curves obtained with Af w = 5 (full triangles), J\f w = 7 
(full circles) and J\f w = 10 (full squares) with periodic bound- 
ary conditions, at Q = 4 ■ 10 — 4 and TV = 50. Open circles are 
MonteCarlo results with Af w = 7. 
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Figure 3. Comparison between the results of the Fourier- 
truncated model and the direct simulation with standard soft- 
ening. The results of the former have been suitably shifted (see 
text). Here g = 4 • 10 -4 and the boundary conditions are reflect- 
ing. Fourier-truncated model with J\f w = 7: TV = 25 (open trian- 
gles), TV = 50 (open circles), TV = 100 (starred polygons). Direct 
TV-body simulation with softening parameter r\ = 0.01: TV = 500 
(starred squares), TV = 1000 (open squares), TV = 2000 (open pen- 
tagons). MonteCarlo canonical results for the Fourier-truncated 
model: TV = 50 (full circles), TV = 100 (full squares). 




E/N 5/3 

Figure 4. Comparison between the results of the Fourier- 
truncated model and the direct simulation with standard soft- 
ening. The results of the former have been suitably shifted (see 
text). Here q = 1 and the boundary conditions are reflecting. 
Fourier-truncated model with J\f w = 7: TV = 25 (open triangles), 
TV = 50 (open circles), TV = 100 (starred polygons). Direct TV- 
body simulation with softening parameter r] = 0.01: TV = 100 
(open squares), TV = 500 (starred squares). MonteCarlo canoni- 
cal results for the Fourier-truncated model: TV = 50 (full circles), 
TV = 100 (full squares). 



the canonical ensemble; however, in the case of a two dimen- 
sional model, it has been observed (Lynden-Bell & Lynden- 
Bell 1977) that the appearance of negative specific heat in 
the microcanonical ensemble corresponds to a signal of a 
phase transition in the canonical ensemble. Nothing similar 
is found here. This is a remarkable result: though a-priori 
the ensemble inequivalence was expected, such a radical loss 
of any signal of the transition in the canonical ensemble was 
unexpected. 



3.3.2 Specific heat 

From the above given definition of temperature the micro- 
canonical specific heat is computed according to the formula 

1 _ dT{E) 
CV ~~ dE ' 

Hence the knowledge of the caloric curve implies the knowl- 
edge also of the specific heat, but in practice the points on 
the T — E plane are inevitably affected by a "noise" which 
can make too rough the computation of CV as the numerical 
derivative of the caloric curve. The interest of an indepen- 
dent numerical derivation of the specific heat is therefore 
twofold: on one side, it constitutes a necessary and useful 
check of the previous results, on the other side it should give 
more information about the order of the phase transition. 
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Figure 5. Fourier-truncated model. Specific heat at N = 50, 
q = 1 and J\f w = 7. Results obtained with reflecting boundary 
conditions (full circles) are compared with those obtained with 
periodic boundary conditions (starred squares). 



Figure 6. Fourier-truncated model. Specific heat obtained with 
(> = 4 ■ 10 — 4 and J\f w = 7. Reflecting boundary conditions. Com- 
parison among: N = 25 (open triangles), N = 50 (open circles) 
and N = 100 (starred polygons). 



The already mentioned Laplace transform technique 
(Pearson et al. 1985) yields the following expression 



cv 



(13) 



iV-2 



which has been used in our numerical simulations; K is the 
total kinetic energy. 

The specific heat vs energy qualitatively displays a pat- 
tern in full agreement with what can be guessed from the 
caloric curves. At small energy values a branch of negative 
specific heat is found. At an energy value slightly below the 
minimum of the caloric curve, the specific heat shows two 
marked spikes, stressing the occurrence of a sudden jump be- 
tween negative and positive values. The high energy asymp- 
totic value is 3/2, in quantitative agreement with the high 
energy slope of 2/3 of the caloric curve. Again the data col- 
lapsing is obtained by scaling the total energy with N 5 ^ 3 . 

In Fig. ^ the specific heat obtained with periodic and 
reflecting boundary conditions respectively is reported vs en- 
ergy density for the Fourier-truncated model. Also through 
this observable it is confirmed that there is a weak sensitiv- 
ity of the outcomes upon the boundary conditions, at least 
as far as thermodynamical observables are concerned. 

For the Fourier-truncated model, Fig. ^ and Fig. |^ show 
- at different values of the density - how the pattern of 
cv changes when N is increased. The corresponding results 
obtained for the direct iV-body simulations are reported in 
Fig. |§] and Fig. ^| A common feature of these results is the 
ambiguous behaviour of the height of the peaks of cv, close 
to the transition point, when N is increased, at variance 
with other systems with short-range interactions (Caiani et 
al. 1998; Cerruti-Sola et al. 2000). The doubtful absence of 
a systematic increase with TV of the peak of cv could be due 
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Figure 7. Fourier-truncated model. Specific heat obtained at g = 
1, Af w = 7 and with reflecting boundary conditions. Open circles 
refer to N = 50 and starred polygons to N = 100. 



to three main reasons: i) a very fine tuning of the energy 
value could be necessary and the right energy values could 
have been missed; ii) the transition could be too mild, for 
example with only a logarithmic divergence of cv\ in) the 
pattern of T(E) could be "V"-shaped, thus d 2 S/dE 2 should 
make a finite jump at the transition point, with the physical 
consequence that the system should undergo a first order 
phase transition. 
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Figure 8. Direct Af-body simulations. Reflecting boundary con- 
ditions and q = 4 • 10 -4 and r\ = 0.01. Comparison of the specific 
heat obtained at N = 500 (starred squares), TV = 1000 (pen- 
tagons) and N = 2000 (starred polygons). 



Figure 10. Comparison of the specific heat obtained at g = 1 
and with reflecting boundary conditions: Fourier- truncated model 
with N = 100 (starred polygons), direct TV-body simulations with 
N = 100 (open triangles) and with N = 500 (open squares). 
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Figure 9. Direct Af-body simulations. Reflecting boundary con- 
ditions, 8=1. Comparison of the specific heat obtained at 
N = 100 and rj = 0.01 (open squares), N = 500 and rj = 0.01 
(open circles), N = 500 and rj = 0.02 (starred squares), and 
N = 2000 and rj = 0.02 (starred polygons). Here rj is varied as a 
stability check on the results. 



The high energy branches of cv, above the transi- 
tion energy, obtained with the Fourier-truncated model and 
with the direct A^-body simulation are in very good agree- 
ment, whereas the quantitative agreement is less good in 
the branch of negative values. Notice that in the Fourier- 
truncated model there is a lower bound to the energy values 
which depends on the number of modes retained in the trun- 



cated Fourier expansion, the larger the number of modes, 
the lower is the minimum attainable value of negative ener- 
gies^. Figures ^ through ^| show that the results of the com- 
putations of cv according to Eq.(^) are in very good agree- 
ment with the corresponding patterns of T(E). In both the 
direct newtonian simulations and in the Fourier-truncated 
model, a kink in the negative slope region of T(E) is found 
and is confirmed by the cv(E/N 5 ^ 3 ) pattern close to the 
transition point. We attribute the appearance of this fea- 
ture to the competing effects of the fixed size of the box 
and the increasing (with energy) gravitational radius of the 
system. 



3.3.3 Order parameter 

Both in the theoretical and numerical studies of phase tran- 
sitions, the choice of a good order parameter for a given 
system is a main point. An order parameter is typically a 
collective variable that bifurcates at the transition tempera- 
ture T c . A classic example is the spontaneous magnetization 
of a ferromagnetic material: it is zero above T c , then it sud- 
denly bifurcates away from zero at T c and increases as tem- 
perature is lowered below T c . True singularities exist only 
in the N — + oo limit; at finite N they are smoothed but it 
is still possible to detect the existence of a phase transition 
by changing iV and observing if the smoothed signals tend 
to sharpen or not. 

In many cases the definition of an order parameter is 



35 Also for a softened newtonian system there is a lower bound on 
the total energy. With the values of the softening chosen here, the 
lower bound is much smaller than that of the Fourier-truncated 
model. 
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not trivial. A-priori this is true also for self-gravitating TV- 
body systems. However, let us notice an interesting property 
of the Hamiltonian (p|) and of the equations of motion (|^): 
at large TV the coefficients S, ^ can be approximated by 
the Si t m, n with an error 0(1/N), and the replacement of the 
Si, m , n by the averages {Si t m,n), computed through a consis- 
tency equation (see Appendix B), decouples all the degrees 
of freedom. The Hamiltonian (^) is approximated by 

.v 

H mf = ^2 + ^ + 

i=l 

16 N ^ m (S ) 

r/ s no , ''T*. — 5T sin(fcia;i) sin(fc m y;) sin(fc n z l ) . 

nL ^ ^ (P + m 2 + n 2 ) 

i—1 Z,m,n— 1 

(14) 

which is now in the form of a sum over independent degrees 
of freedom, because the (S;, m ,„) are collective variables or, 
equivalently, order parameters, that take the same values 
for all the coordinates Xi,yi,Zi. The replacement of the co- 
efficients 5 , (, mi „ with the averages (5 , (, mi „) in Hamiltonian 
( [li] ) is a typical mean-field approximation in statistical me- 
chanics (the (S;, m ,n) are like Weiss molecular field of the 
early times of statistical mechanics of magnetic materials). 
The recasting of Hamiltonian ([!]) into the approximate form 
(p^), which can be as precise as we may want if TV, Af w 
and L are sufficiently large, is of great prospective interest 
for theoretical - analytic or semi-analytic - computations of 
statistical mechanical kind (see Appendix B). 

Out of the huge family of order parameters (Si tm , n ), it is 
natural to define the following two global order parameters 
in analogy with other more standard contexts (Caiani et al. 
1998; Cerruti-Sola et al. 2000) 

M 1 = i V + % (g,, m ,n) (15) 

TV ^ (P+m 2 +n 2 ) 

l,7n,n — l 

and 

M 2 = i £ (P + £ + n») ■ (16) 

l,m,n — 1 

The coefficients (? + m + n)/(/ 2 +m 2 +ra 2 ) are the sum of the 
corresponding coefficients in the equations of motion ^ for 
the three variables x, y, z. They decrease as the norm of the 
wavevector k — (I, m, n) increases, compensating the grow- 
ing number of the Si >m , n terms as the norm of k = (I, m, n) 
increases. We easily realize that Mi cannot be very sensitive 
to a change of the particle distribution in the box because 
the dishomogeneities measured at different scales tend to 
cancel each other. Thus, to get rid of this problem, M2 is 
defined through the sum of the absolute values of the S;,m,n- 
Actually, the numerical computations, performed at differ- 
ent energies and TV, confirmed that Mi is always very close 
to its minimum value (which depends on Af w )- At variance, 
the order parameter M2 is very effective to detect the clus- 
tering transition and to make strict the analogy between this 
transition and a thermodynamic phase transition. 

In order to measure the deviation from the homoge- 
neous spatial distribution of particles, it is convenient to 
slightly modify the order parameter as follows: we denote 
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Figure 11. Order parameter /J2 for the Fourier-truncated model. 
g = 1, Afw = 7 and reflecting boundary conditions. Open circles 
refer to TV = 50, starred polygons to TV = 100, starred squares to 
TV = 200 and open squares to TV = 400. 



with Mo the value of M2 that corresponds to a uniform 
density of particles in the box, then we modify the order 
parameter M2 to 112 = v(Mi — Mo), where v is a normal- 
ization constant. Thus /X2 varies in the interval [0,1], be- 
ing zero for a uniform distribution of particles and one for 
a maximally clustered configuration. Such a normalization 
makes easy the comparison of the results obtained at dif- 
ferent TV. In the case of the Fourier truncated model, the 
maximally clustered configuration is obtained when all the 
particles are concentrated in a cell whose side is equal to the 
shortest wavelength. However, for simplicity, we empirically 
normalized ^2 by using the largest value of M2 measured at 
the lowest accessible energy for the given truncation of the 
Fourier expansion. 

The values Mo of A/2, corresponding to a uniform den- 
sity of particles in the box, have been numerically com- 
puted by averaging Mi on a few hundreds of uniformly 
distributed random configurations of 50, 100, 200 and 400 
particles. For TV = 50, 100, 200, 400 and Af w = 7 we found 
M = 3.49, 2.60, 1.97, 1.54 respectively. 

The averages (\Si t m,n\) are numerically computed as 
time averages to give M2, and the order parameter /12 is 
worked out as a function of the energy and of the number 
of particles by using the above reported values of Mo. 

The same order parameter jj,2 has been computed for 
the direct TV-body simulations. Also in this case we consid- 
ered Af w — 7, though in principle this choice is arbitrary 
because Afw is independent of the dynamics. Since a phase 
transition is a collective effect, driven by the long wavelength 
modes, there is no need for a large value of Afw ■ 

The final results are reported versus E/N 5/3 in Fig. 
for the Fourier-truncated model, and in Fig. [l^ for the direct 
TV-body simulations. 

The data reported in these figures display typical pat- 
terns that are found in presence of a phase transition in other 
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12. Order parameter p,2 for the direct TV-body simula- 
= 4 • 10 — 4 , Af w = 7 and reflecting boundary conditions. 
= 2000. 



contexts. In fact, the order parameter fj,2 displays a classic 
bifurcation pattern at the transition point, smoothed by the 
finite and not very large number of particles but with a clear 
tendency to become sharper at larger N as is shown by Fig. 
jljj . The simulation of the direct iV-body system allows to 
work out the bifurcation pattern of the order parameter in 
a larger energy interval, because there is no lower bound to 
the energyP] The inset of Fig. [l^ displays the pattern of (12 
on a large energy interval. The vertical dotted lines in both 
figures indicate the energies of the peaks of specific heat. In 
general, the inflection point of the order parameter and the 
location of the peak of the specific heat do not coincide but 
tend to get closer at larger N, thus making more precise the 
phase transition point. The bifurcation pattern of \i2 gives a 
reliable indication about the order of the phase transition: it 
strongly indicates a second order transition, thus eliminat- 
ing the ambiguity of the information given by cv. In fact, 
in presence of a first order transition, the order parameter 
/i2 should make a finite jump at the transition point which 
is clearly not our case. We have formulated in Section 3.3.2 
two possible explanations (i) and iij) - alternative to the 
first order phase transition - of the behaviour of cv which 
can reconcile the outcomes on cv with those on /X2. 

Notice that M2 is also a measure of the average modulus 
(\Fi\) of the force acting on each mass and - independently 
of N - this is never negligible for a bounded and clustered 
system, whereas in the homogeneous phase depends 
on the relative density fluctuations and thus it is expected 
to decrease with increasing N. The physical consequence is 
that, in the homogeneous phase, the larger N the smaller 
is the relative average weight of the potential energy with 



^ As already noticed, the smoothing parameter entails a lower 
bound to the energy also in this case, but it occurs at a very large 
negative value. 



respect to the kinetic energy so that the system behaves like 
a collection of almost independent particles, i.e. not far from 
a perfect gas. 



3.4 Equation of state 

The equation of state of a system contains the basic infor- 
mation of how pressure changes with volume at constant 
temperature. This information is contained in a family of 
isothermal curves in the P—V plane. If a phase transition is 
there, then the isotherms at T < T c display a substantial dif- 
ference from those at T > T c , as an example, the isotherms 
of the van der Waals equation, describing the liquid-vapour 
transition, display a kink at T < T c (Huang 1963). 

Since dynamics has its natural statistical counterpart 
in the microcanonical ensemble, where energy rather than 
temperature is the independent variable, the isotherms are 
replaced by isoenergetic curves in the P — V plane. 

From thermodynamics we have 



/dS\ _ P 
\dVj E ~ T 



and using the above given definition of temperature one gets 
fcflJ = , 

where f2jv is defined in the Appendix A. Hence 
8Q.N * 



P = 



( 1 dn N \ 

\ujn dV J j 



whence, with the aid of the already mentioned Laplace trans- 
form technique, one obtains (Pearson et al. 1985) 



= y(k B T) u 



(17) 



which is the microcanonical equation of state yielding pres- 
sure vs volume at constant energy; K is the total kinetic 
energy and U is the potential energy. Temperature is given 
by EqjlfJ) as a function of the energy. Again, the averages in 
Eq. (117) are computed as time averages along the numerical 
phase space trajectories. 

By replacing L = V 1 / 3 in the potential part of the 
Hamiltonian (|^), and noting that after the replacements 
Xi = Lxi, Vi = Lrji and z t = LQ (where %i, are adimen- 
sional variables) the arguments of the sines do not depend 
on L, the following simple result is found 



dU _ _U_ 
W ~ ~3V ' 

so that the microcanonical equation of state finally reads 



(18) 



N 

T 



3V 



{K-X 



(19) 



Its explicit analytic computation seems hard but still feasi- 
ble with the aid of the mean-field decoupling of the degrees 
of freedom in Eq. (|l4|) , however such an attempt is beyond 
our present aims and again we use dynamics in order to es- 
timate the microcanonical averages through time averages. 
Moreover, by extracting out of U its prefactor V"~ 1//3 [see 
Eq.(§)], we see that Eq.([l9|) is of the form 



PV = Nk B T -U(E)V 



-1/3 



(20) 
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where 14(E) would be the outcome of the analytic com- 
putation of the microcanonical averages, i.e. U(E) = 
l(V 1/3 U/K)(K' 1 }. Equation agrees with the already 
proposed modification of Boyle's perfect gas law for self- 
gravitating systems (see the next Section) which is here qual- 
itatively rederived in the microcanonical ensemble and quan- 
titatively computed numerically. The results are reported in 
Fig.jl3|, where we have differently marked the points that cor- 
respond to positive and negative specific heat respectively. 
The one-to-one correspondence polytrope-negative specific 
heat (perfect gas law-positive specific heat) is well evident. 
The transition from a polytropic V -4 / 3 law to a V~ x perfect 
gas law is well evident. This reflects the competition between 
the two terms in the r.h.s. of Eqs. (jT|) and (§o|). Figure [lj 
shows PV/N 2 vs V, obtained at N = 100, 200 and at dif- 
ferent values of the energy E which have not been shifted 
towards their corresponding newtonian values as it was done 
in the preceding Sections (the reason should be clarified by 
the concluding remarks of the present Section). In order to 
compare the results obtained at different N, E is scaled with 
iV 2 as well as PV . The results of Fig. ^ directly compare 
to Eq.(pc|) and make clearer the cross-over between V -4 / 3 
and V~ . Such a cross-over seems rather smooth, though a 
discontinuity in the derivative dP/dV cannot be excluded, 
and this seems to confirm that the gravitational clustering 
phenomenon is a peculiar phase transition with respect to 
all the known laboratory phase transitions. 

The non-trivial physics behind the cross-over is phe- 
nomenologically displayed by Figs. [IB] through [^j, where a 
few snapshots of the spatial distribution of the N interacting 
masses are projected onto the x — y plane. It turns out that 
the I/ -4 / 3 branch of the equation of state corresponds to the 
clustered phase, whereas the same picture obtained in the 
V -1 branch displays an homogeneous distribution of parti- 
cles, typical of a gas phase. Figures |l5| and [l^ refer to the 
Fourier-truncated model and are obtained keeping energy 
constant and varying the volume of the box, whereas Figs. 
[It] through [llj refer to the direct JV-body system and have 
been obtained by keeping the volume constant and varying 
the energy. 

A remarkable property of this equation of state is the 
absence of a critical point of flat tangency, i.e. dP/dV = 0, as 
it occurs in the liquid-gas transition described by the van der 
Waals equation. The kink in the van der Waals isotherms 
below the critical one (corresponding to the phase coexis- 
tence and to the existence of the latent heat) is here absent, 
thus the "gravitational condensation" appears rather differ- 
ent from the gas-liquid condensation. 

A final comment about the energy values of the P — V 
curves is in order. At any finite JV W , the Hamiltonian (^]) 
implies a truncation error in the potential energy with re- 
spect to the Hamiltonian ([l]). At fixed volume this implies 
just a shift in energy scale, making easy the comparison be- 
tween the Fourier-truncated model and the direct model. 
But when volume is varied, also the truncation error AU 
varies. It is numerically confirmed that the expected In- 
dependence of the correction to the potential energy scale is 
proportional to V _1//3 (for example, at N = 100, A/" = 7, 
V = 10 2 and V = 2.5 10 s , we found AU/N 5/3 = -0.85 
and AU/N 5/3 = -0.065 respectively). Thus, the trunca- 
tion can only shift the value of V at which the crossover 
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Figure 13. Equation of state: pressure vs. volume. Fourier- 
truncated model with N = 100, Af w = 7 and reflecting boundary 
conditions. Open triangles refer to E = 3174.8, full and open cir- 
cles to E = 793.7. E are not shifted (see text). Full circles denote 
points that correspond to negative specific heat, conversely open 
circles denote positive specific heat. Reference lines (long dashed) 
correspond to a polytropic V~ 4 / 3 law and to a perfect gas V - 1 
law respectively. 



occurs, without qualitatively changing the phenomenology. 
The larger A/" the smaller the shift will be. 



4 COMPARISON WITH THEORETICAL 
PREDICTIONS 

It was suggested by Terletsky (1952) that for a large mass 
M of gas in a volume V and at temperature T, composed 
of N particles under a boundary pressure p, the equation of 
state of the perfect gas (Boyle's law) should be corrected to 



PV = Nk B T - aGM'V 



-1/3 



(21) 



where a accounts for the shape of the mass and G is the 
gravitational coupling constant. Shortly after, it was shown 
that such a correction of Boyle's law results in gravitational 
instability (Bonnor 1956). 

Two theoretical milestones followed with the papers of 
Antonov (1962) and Lynden-Bell & Wood (1968), where 
a gravitating system of point particles in a spherical box 
and the thermodynamics of a self-gravitating isothermal 
gas sphere were respectively considered. Antonov's very in- 
teresting result is about the existence of the gravothermal 
catastrophe when the density contrast g c / g e between centre 
and edge of the box exceeds the value 709. Lynden-Bell and 
Wood (1968) related Antonov's result with the existence of 
negative heat capacity and found the remarkable result that 
stable isothermal spheres with negative specific heat exist in 
the density contrast range between 32.2 and 709; the specific 
heat diverges at Qc/Qe = 32.2 and becomes positive below 
this threshold value; above g c / g e — 709 there is a runaway in 
the catastrophic phase. Moreover, the minimum attainable 
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Figure 14. Equation of state for the Fourier-truncated model 
with TV = 100, Mw = 7 and reflecting boundary conditions. 
PV/N 2 is plotted vs. V to put in evidence the transition from 
a polvtropic law to a perfect gas law. Symbols are the same as in 
Fig. |l§ Starred polygons refer to TV = 200 and E = 3174.8. 



Figure 16. Superposed snapshots of particle positions projected 
onto the x — y plane for the Fourier-truncated model with TV = 
100, Af w = 7 and reflecting boundary conditions. Energy E = 
793.7 and log V = 2.0. Clustered phase (see Fig. |l4| ). 





Figure 15. Superposed snapshots of particle positions projected 
onto the x — y plane for the Fourier-truncated model with TV = 
100, Af w = 7 and reflecting boundary conditions. Energy E = 
793.7 and log V = 5.0. Homogeneous phase (see Fig. 113). 



Figure 17. Snapshot of particle positions projected onto the x — y 
plane for the direct TV-body system with TV = 2000 and reflecting 
boundary conditions, g = 4 ■ 10 -4 , T3/TV 5 / 3 = —0.5 and logV = 
6.7. Clustered phase (see FigJ3J). 



temperature by an equilibrium isothermal sphere of radius 
r e and mass M is 



GmM 
2.52k B r e 



(22) 



and is achieved for a density contrast Q c / Qe = 32.2 where 
cv becomes negative. 



Our present study naturally represents the microscopic 
counterpart of the thermodynamic framework where all the 
above mentioned results were obtained. 

In the preceding section we have seen that the 
dynamical-microcanonical equation of state of a self- 
gravitating TV-body system yields a numerical result in per- 
fect agreement with that of Eq.(|2l]), at least as far as the 
functional form of the equation of state is concerned, and 
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Figure 18. Same as Fig.^but with E/N 5 / 3 = -0.03. Transition 
regime (see Fig.H). 



Figure 20. Density contrast vs energy per particle for the 
Fourier-truncated model with N = 100, Af w = 7 and reflecting 
boundary conditions. The vertical dashed line separates cy < 
(left) from cy > (right). The horizontal dashed line indicates 
the density contrast centre to edge Qc/Qe = 32. 



50 - 



>> - 



-50 - 




Figure 19. Same as Fig.^ but with E/N 5 / 3 = 0.12. Homoge- 
neous phase (see Fig.H). 



even though we obtained isoenergetic P — V curves instead 
of isothermal ones as is implicit in Eq.(^). The transition 
point, separating Boyle's law from the polytropic one, actu- 
ally corresponds to the loss of stability of the homogeneous 
gas phase which goes into a clustered state, as is also quali- 
tatively shown by Fig.s ^| and [H| Moreover, we have found 
that the transition actually occurs at the minimum attain- 
able temperature^ in fairly good agreement (taking into 



III We have G = m = kg = 1, r e = L and M = N. Trivial al- 
gebra gives T mln /N 2 / 3 = (N/V) 1 / 3 /2.52, whence T mln /N 2 / 3 ~ 



account the difference of the geometry of the box) with the 
prediction of Eq. (^) , and that below the energy that corre- 
sponds to the minimum attainable te mpera ture the specific 
heat is negative, as shown in Section 3.3.2. The increase of 



the density N/V implies a shift of the transition energy den- 
sity towards more negative values - as shown by Figs. |§] and 
^ - and this is coherent with the Lynden-Bell & Wood pre- 
diction of a critical density contrast above which cy < and 
clustering occurs. In fact, increasing the density amounts to 
decrease the density contrast (if E/N is kept constant) and 
increasing E/N amounts to decrease the density contrast 
(at constant q), thus, in order to keep constant a given ratio 
Qc/Qe (for example the threshold value) it is necessary to 
decrease E/N. 

Finally, in Fig.^, we show that the density contrast 
between centre and edge when cy becomes negative is ac- 
tually in strikingly good agreement with the prediction of 
Lynden-Bell & Wood. Even if these authors warned that 
"both these instabilitie s]* * *| depend on the presence of the 
heat-bath; thermally isolated systems do not suffer this type 
of instability" , we found that the clustering instability shows 
up also in the absence of a heat-bath, whereas this is not the 
case of Antonov's gravothermal instability. 

There are several possible physical explanations for the 
absence of a dynamical counterpart of the gravothermal 
catastrophe. A first possibility is that its growth rate is 
so small that a huge time is required to observe it; this is 
what has been found for a hybrid model - mean-field with 



0.4 for Q = 1, whereas we found 0.6 - as shown in Figure [ij — and 
Tmin/N 2 /' 3 ~ 0.03 for g = 4 10~ 4 , whereas we found 0.04, as 
shown in Figure |i| 

* * * Occurring at Qc/Qe = 32.2 and at the Antonov threshold 
Qc/Qe = 709. 
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local perturbations - where the time scale to develop the 
gravothermal catastrophe is larger than the Hubble time 
(Sygnet et al. 1984). Another dynamical mechanism that 
can inhibit the collapse is the formation of binary systems 
(Heggie & Aarseth 1992). 



5 CONCLUSIONS 

The present paper has mainly dealt with the dynamic ( "mi- 
croscopic") origin of the thermodynamic behaviour of self- 
gravitating systems and - in close connection with this topic 
- with the problem of their appropriate statistical mechan- 
ical treatment. 

The dynamics of self-gravitating iV-body systems, de- 
scribed by two different kinds of regularization of the newto- 
nian interaction, has been numerically investigated. Dynam- 
ics provides the time averages of certain quantities which - 
with the aid of suitable formulae developed in the field of 
Molecular Dynamics - yield thermodynamics. The link be- 
tween dynamics and thermodynamics is made in the micro- 
canonical ensemble of statistical mechanics. 

The regularizations of the newtonian interaction are the 
standard softening and a truncation of the Fourier expansion 
series of the two-body potential. N particles of equal masses 
are constrained in a finite three dimensional volume. 

The introduction of a Fourier-truncated model makes 
possible a mean-field-like decoupling of the degrees of free- 
dom, which is of great prospective theoretical interest for a 
statistical mechanical treatment of self-gravitating systems. 
Moreover, just through such a decoupling, an order parame- 
ter can be naturally introduced to signal the occurrence of a 
phase transition. Actually, through the computation of the 
caloric curve, of the specific heat and of the order parameter, 
clear evidence is found for a dynamical signature of a phase 
transition of clustering type which appears as a second or- 
der phase transition. Thus, the gravitational condensation 
seems to take place trough a transition milder than that of 
a gas- liquid condensation, which is a first order transition 
with a finite jump in the entropy. This is also coherent with 
our numerical computation of the microcanonical equation 
of state: the counterpart of the clustering transition on the 
P — V plane is a cross-over between a polytrope of index 
n = 3, i.e. P = KV~ 4/3 and a perfect gas law P = K'V' 1 , 
without any kink in the P vs V curves as it would be ex- 
pected for a first order transition. 

The very good agreement between direct simulations of 
the TV-body systems and Fourier-truncated models is also 
interesting because it reveals that the relevant informations 
for thermodynamics are mainly conveyed by the large spatial 
scales rather than by the small ones. This is in agreement 
with the fact that the statistical mechanical peculiarities of 
TV-body self-gravitating systems are due to the long-range 
nature of the newtonian interaction and are not influenced 
by its short-scale singularity. 

The scale invariance (Heggie & Mathieu 1986) of the 
gravitational TV-body systems is broken by constraining the 
particles in a box. The physical reason is that the box in- 
troduces a new length scale besides the gravitational one. 
The clustering phase transition is a consequence of the ex- 
istence of these two independent length scales. Therefore, 
the physical meaning of what we found is that whenever an 



extra length scale is present besides the gravitational one, 
an TV-body system can undergo a clustering transition. The 
use of the box is perhaps the simplest way of breaking scale 
invariance, but a halo of diffuse matter or any other external 
potential could in principle work. 

We have found that newtonian dynamics naturally 
yields a number of classical results obtained within a phe- 
nomenological (thermodynamic) framework. A remarkable 
difference exists between the theoretical prediction of two 
kinds of transitions, the clustering and the gravothermal 
catastrophe, and our numerical results: only clustering is re- 
covered in the dynamical approach. Actually, the possibility 
of a dynamical suppression of the gravothermal catastrophe 
(motivated either by a huge instability time or by the for- 
mation of binaries) has been discussed in the more recent 
literature. 

Finally, we have compared the microcanonical (dynam- 
ical) averages to their canonical ensemble counterparts, ob- 
tained through standard Monte Carlo computations. A re- 
markable result found here concerns the ensemble inequiva- 
lence. This is a-priori expected because of the appearance of 
negative specific heat - forbidden in the canonical ensemble 
- in the dynamically worked out thermodynamics. However, 
the inequivalence is so strong that no signal of the cluster- 
ing transition seems to survive in the canonical ensemble. 
In the case of other models some trace is left (Lynden-Bell 
& Lynden-Bell 1977), but perhaps the transition here is too 
soft. A reliable statistical mechanical approach to (regular- 
ized) self-gravitating systems seems possible only in the mi- 
crocanonical ensemble. 
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APPENDIX A: 

For the sake of clarity and self-containedness, in this Section 
we briefly recall some basic facts about dynamics, statistical 
mechanics and thermodynamics. 
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Al From dynamics to statistics 

The process of measuring an observable f(q,p), function of 
the microscopic coordinates q and momenta p of the parti- 
cles composing a macroscopic system, always conceptually 
involves a time average, necessarily carried out in a finite 
time interval T, of the form 



h{t) 



t+T 



f(q{r),p(T))dT 



(Al) 



Let us now suppose that the system is at equilibrium, i.e., 
that the time average (Al), for sufficiently large T, does not 



depend upon t, whence we can write 



/= lim f T (t) 

T— >oo 



(A2) 



The fou ndi ng hypothesis of equilibrium statistical mechanics 
is that ( A2 ) exists and is independent of the initial condition 
(q(to) , p(to)) , except for a zero measure set in the accessi- 
ble region of phase space. In the case of an isolated system, 
such a region is the constant-energy hypersurface Eb, and 
the measure we are referring to is the dynamically invari- 
ant measure in the phase space — the Liouville measure — 
restricted to E_g. If this hypothesis, known as the ergodic 
hypothesis, is valid, then for any observable f(q,p) we have 

7=(/>m, (A3) 
where (-) M denotes the microcanonical phase average 



(/>. 



dqdpf(q,p) 5[H(q,p)-E} : 



(A4) 



where p — {pi, . . . ,pn} and q = {qi, . . . ,gjv} and lu is the 
normalization 



ojn(E) 



dq dp S[H(q,p) - E] 



(A5) 



The erg odic hypothesis makes the problem of the com- 
putation of ( |A2j ) no longer dependent on the knowledge of 
the dynamics. If we consider a non-isolated system, in con- 
tact with a much larger system, and we apply Eq. (AT) to 



the sum of the two systems, we obtain the canonical formu- 
lation of statistical mechanics. As is well known, the Gibbs' 
probability density in phase space is now proportional to 
cxp(-f3H) rather than to S(H - E), where f3 = l/k B T (k B 
is the Boltzmann constant and T is the temperature), so 
that 



</>c 



where 



1 

Zn 



dp f(q,p) exp[-(3H(q,p)] , 



Z N {T,V)= / dq dp exp[-/3H(q,p)] 



(A6) 



(A7) 



is the canonical partition function. 

It is natural to think of dynamics as the basic source 
of statistical mechanics, and therefore one can wonder un- 
der what conditions the dynamics is ergodic, so that Eq. 



( A5 ) holds. Even if there is no rigorous proof of ergodicity 



* * * Szilard (1925) showed that the functional form of the canon- 
ical distribution is constrained by a consistency requirement with 
the Second Law of Thermodynamics. 
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but for abstract systems, there are only trivial examples of 
physical systems where it is possible to observe a violation 
of ergodicity (with the exception of amorphous materials, 
like glasses or spin-glasses), i.e. where the predictions of sta- 
tistical mechanics are in disagreement with the outcomes of 
experimental measurements at equilibrium. 

The real problem is the dynamical realization of equi- 
librium. In statistical mechanics it is ab initio assumed that 
a system is at equilibrium, but the equilibrium concept is 
very subtle, and depends in a crucial way on the time scales 
over which a phenomenon is observed. In Feynman's words, 
at equilibrium all the "fast" things have already happened, 
while the "slow" ones not yet. 

This is the reason why from a physical point of view it is 
mixing, rather than ergodicity, the most important feature of 
a dynamical system. A mixing dynamics is also ergodic; the 
converse is not true. Under a mixing dynamics, a generic 
probability distribution in phase space o(q,p,t) converges 
exponentially fast to the stationary equilibrium distribution. 
Therefore mixing is responsible for the relaxation to equilib- 
rium, and for the convergence of time averages to statistical 
equilibrium values on finite time scales. This sheds a new 
light on the role of dynamics in explaining thermodynami- 
cal phenomena, because there is a tight relationship between 
mixing and that kind of dynamical instability which is called 
deterministic chaos. 

The reason to believe that chaos is at the origin of phase 
space mixing is that in all the systems where the mixing 
property can be rigorously ascertained, mixing is found to 
be a consequence of chaos. 

A2 From statistics to thermodynamics 

Statistical mechanics also provides a link between the mi- 
croscopic (dynamic) and macroscopic (thermodynamic) de- 
scriptions of large collections of objects. The main frame- 
works are (Huang 1963) the microcanonical ensemble (when 
energy and number of particles are given), the canonical 
ensemble (when temperature and number of particles are 
given) and grand-canonical ensemble (when total energy and 
number of particles are allowed to fluctuate). The equiva- 
lence of these ensembles, in the large TV limit (thermody- 
namic limit), is a fundamental point, so that the choice of 
the statistical ensemble is only a matter of practical conve- 
nience. 

In the microcanonical ensemble the link with thermody- 
namics is made through the following definition of entropy 
(V is the physical volume, TV is the number of particles) 



S(E, V,N) = k B log lo n (E). 



(A8) 



Another definition of the microcanonical ensemble volume 



in phase space, alternative to u>n in Eq. (A5), is 



Sl N (E)= / dq dpO[E- H{q,p)] 



(A9) 

where 0(-) is the Heaviside step function, and the entropy 
is now given by 

S(E,V,N) = k B log Q N (E). (A10) 

The two definitions of entropy differ by 0(1/N) terms and 
therefore coincide in the thermodynamic limit. With both 
definitions of entropy, temperature is obtained as 



T(E) 



dS_ 

dE 



(All) 



Strictly speaking, the above definitions of entropy are 
arbitrary and justified a-posteriori mainly by showing that 
they are consistent with the laws of thermodynamics (Huang 
1963). 

In the canonical ensemble the link with thermodynam- 
ics is based on the following definition of the Helmoltz free 
energy 



F{T,V,N) = -k B T log Z N (T,V) , 



(A12) 



where Zn is the partition function defined in Eq. ( A7); from 



the function F all the other thermodynamic functions are 
obtained through standard Maxwell's relations. Also the def- 
inition of F(T, V, TV) is a-priori arbitrary and is given validity 
a-posteriori. 

At the macroscopic level, a many-body system has a 
good thermodynamic behaviour only if the microscopic in- 
teraction potential fulfils two basic properties: stability and 
temperedness (Ruelle 1969). Loosely speaking, temperedness 
means that the positive part of the interaction energy be- 
tween particles at large distances is vanishingly small. The 
negative part of the interaction is left arbitrary by the def- 
inition of temperedness, thus it has to be accompanied by 
the stability condition - which is the relevant condition for 
gravitational interactions - 



,rjv) > —NB 



(A13) 



where B > 0. From these conditions the existence of the 
thermodynamic limit and the equivalence of the statistical 
ensembles follow. The existence of the thermodynamic limit 
means that by letting TV — > oo and V — > oo so that N/V 
remains finite, then also the energy density, entropy den- 
sity and Helmoltz free energy density remain finite. In other 
words, energy, entropy, free energy are extensive quantities. 

The ensemble equivalence holds in the thermodynamic 
limit and means that the same thermodynamics can be de- 
scribed by means of microcanonical and canonical ensem- 
bles, for example. 



APPENDIX B: 

Let us briefly sketch here how the computation of the coef- 
ficients of the mean-field Hamiltonian (jl4|) should proceed 
and, similarly, how this mean-field version of the gravita- 
tional TV-body Hamiltonian (g) would make feasible some 
analytic computation of microcanonical statistical averages. 
The microcanonical average {Si trn , n } of any coefficient 



= y sm(kiXj)sm(k m yj)sm(k n Zj) 



is in principle computable according to Eq.(A4), where 
H(p, q) is given by Eq.(^|). However, in practice such a com- 
putation is unfeasible. If the Hamiltonian (^) is replaced 
with the mean-field Hamiltonian (|l4|), then 



1 f N 

/ TT^i dpi 5[H„ 

UN J AA 



,f{Pk,Xk, (Sl, m ,n)) ~ F] Si 
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This integral equation in the unknown (Si, m . n ) is the con- 
sistency equation for the mean-field approximation. 
Now the following identity is very useful 



c- x U 



1 f N 

/ I I dSi dp S[H m f(p k ,x k , (Si, m ,„)) - E]Si, r , 

i^jv J t 



where C stands for Laplace transform and £ _1 for its in- 
verse. In fact, the momenta are easily integrated (Pearson 
et al. 1985) leaving 



1 LOnCq \sJ 



n 



dxi e aU i@*>( s t,m,n>) sin(fc;j:i) sm(k m yi) sin(k n Zi) 



where Co is a numerical constant, s is the variable of the 
Laplace transform, and thanks to the decoupling of the de- 
grees of freedom in H m f, all the triple integrals in square 
brackets are equal and independent one from another be- 
cause the potential energy U in Eq.(^) is the sum of in- 
dependent contributions Ui. If we denote the generic triple 
integral with F(s, (Si, m , n )), we are finally left with 



which, in general, can be worked out through some approx- 
imate method, the most popular being the saddle point 
method (Huang 1963). 

Once the coefficients (Si. m ,n) have been computed, the 
mean-field Hamiltonian H m f is completely specified and can 
be used to compute other microcanonical averages through 
the same guideline depicted above. 



This paper has been produced using the Royal Astronomical 
Society/Blackwell Science L^TfrjX style file. 
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